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I. INTRODUCTION 

The ability of controlling parts or subsystems of a quantum system is the main resource 
of a future quantum computer with its processing superiority over nowadays computers. 
Such ability presupposes the control of local operations on parts of the system and its 
effects due to quantum correlations on other parts of the system. Such quantum correla¬ 
tions among distinct parts of the system are encoded by the keyword entanglement - the 
characteristic trait of quantum mechanics, in the words of Schrodinger. It is therefore of 
great interest to be able to quantify or measure how “entangled” the parts of the system 
are. This is the reason of the great effort being carried in the past decades to elucidate 
the notion of quantum entanglement. 

Entanglement has been studied in diverse systems, the simplest and most studied 
being those associated with quadratic spinless fermionic chains. These are related via a 
(non-local) Jordan-Wigner transformation to simple spin systems such as Ising, XX and 
XY with an external magnetic coupling. It is fair to say that nowadays we have a good 
control of almost all aspects of quantum correlations in such systems. Such achievement 
were possible due to the sophisticate resolution of diverse mathematical problems, most 
of them associated with the computation of the determinant of correlation matrices - 
usually of the block Toeplitz type due to translational symmetry. The sophistication may 
be appreciated by noting that such computations are instances of the Riemann-Hilbert 
problem (RHP). This fact brings in a geometrical flavour to the study of entanglement 
entropy for such systems. It is a question of the analysis of the topology of a certain 
Riemann surface associated to the system. 

In the framework of this successful ideology, we can find two classes of block Toeplitz 
matrices [1]. On one hand, those which can be reduced to a pure Toeplitz one, with 
scalar symbol, where the Fisher-Hartwig theorem [2J is the key to obtain an asymptotic 
expansion of our determinant. That is the case of spinless fermionic chains described by 
a Hamiltonian which preserves the fermionic number. In this way, Jin and Korepin [3] 
obtained a formula for the von Neumann entropy associated to a set of contiguous sites in 
the XX spin chain. Using the same philosophy, the Renyi entropy of an interval for every 
stationary state of a general quadratic spinless fermionic Hamiltonian which preserves 
the number of excitations is obtained in [JJ. On the other hand, we have those block 
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Toeplitz matrices which cannot be reduced to a scalar one; as usually happens when the 
above number symmetry is broken. In general, we cannot bypass the RHP. Confronting 
it, Its, Jin and Korepin |5] obtained a formula for the von Neumann entanglement entropy 
associated to a set of contiguous sites in the XY spin chain. An expression for the Renyi 
entropy was then written down by Franchini, Its and Korepin [6]. Later the von Neumann 
entropy for the case of quadratic spinless fermionic chains with finite range coupling was 
provided by Its, Mezzadri and Mo [7j. In a different spirit, using a duality transformation, 
Peschcl [Bj obtained the same kind of formula for the XY spin chain. 

Later on, Kadar and Zimboras [2] started a systematic analysis of the entanglement 
entropy for self-dual models with broken reflection symmetry. It is interesting to note 
that in this self-dual case the block Toeplitz matrix analysis reduces to that of a scalar 
Toeplitz matrix. 

In this work, we extend further the above sequence of complexities by considering the 
Renyi entanglement entropy for quadratic spinless fermionic chains with complex finite 
range couplings. Due to the complex nature of the couplings, some new possibilities, to 
our knowledge not known in the literature, appear. We have to deal with some interesting 
challenges associated with some discontinuities in the symbol of a block Toeplitz matrix. 
As a resolution of these challenges we here propose a new formula for the determinant of 


such Toeplitz matrix with 2x2 symbol (23): 


log D x (X) = ® f logdet [XI - g(9)}d9 + log \X\ + • • • , 

11 J ~' K 3=1 

where bj encapsulates the information due to the discontinuities of the symbol. The 
important fact of this formula is its logarithmic scaling with the size of the subsystem 
|Y|. As an application of this new formula, we proceed to an exhaustive analysis of the 
entropy computation of a spin chain with Dzyaloshinski-Moriya (DM) coupling and a 
Kitaev fermionic chain with long range pairing. 

As stated by Its, Mezzadri and Mo in [7], At the core of our derivation of the 
entropy of entanglement is the computation of determinants of Toeplitz matrices for a 
wide class of 2 x 2 matrix symbols.”. The 2x2 matrix character of these symbols stems 
from the non-invariance of the fermionic number. This is a distinctive feature of the 
XY model. We also consider finite range couplings; being the Kitaev chain with long 
range pairing an example. The novelty in our work lies in the complex nature of the 
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couplings. This feature allows one to consider systems with broken reflection and/or 
charge conjugation. An example of a system with broken reflection symmetry is the DM 
spin chain. 

When considering such systems with broken reflection and/or charge conjugation, we 
obtain a wider moduli of system, among those one may encounter systems with critical 
behaviour. Such systems are expected to display an entropy that scales logarithmically 
with the size of the system. This is indeed what we obtain when we use our proposed 
formula for the determinant of the block Toeplitz matrix in the Renyi entropy formula. 

We organize this paper as follows. In section 2, we present the general Hamiltonian for 
a spinless fermionic chain with finite range coupling and define the Fock space of states; 
in section 3, we review the well known relation between the entropy and the correlation 


matrix and in (12) we derive the latter for a general stationary state; in section 4, we 
compute the determinant of the block Toeplitz matrix, first in the smooth case where the 


Its, Mezzadri and Mo’s formula (19) applies, and then by extending it in (23) to the case 
where discontinuities are present; later we use this result to obtain our general expression 


(26) for the asymptotic behaviour of the entanglement entropy in a critical theory with 
long range couplings and the above broken symmetries; in sections 5 and 6, we apply the 
formulae of previous section for an analysis of the entropy in the DM spin chain model 
and in a Kitaev chain with long range pairing respectively; finally, section 7 contains our 
conclusions. 


II. FINITE RANGE HAMILTONIAN 


Let us consider an IV-site chain of size i with N spinless fermions described by a 
quadratic, translational invariant Hamiltonian with finite range couplings (L < N/ 2) 

1 N L 

H ~ 2 EE ( 2 Aia) n a n +i + Bia) n ci n _ ^ R;a n a n _)_^ . (1) 

Z n =1 l=-L 

Here a n and a 1 represent the annihilation and creation operators at the site n. The only 
non-vanishing anticommutation relations are 


o m } 6 n 


and in (|Tj) we assume a n+ N = a n . 
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In order that H is Hermitian, the hopping couplings must satisfy A_i = Ai. In addition, 
without loss of generality, we may take Bi = —B_i. When Bi = 0, VZ, our Hamiltonian 
preserves the fermionic number. If Ai is not real the reflection symmetry (n i—>■ N — n) is 
broken. The imaginary part of Bi breaks charge conjugation (a n h->■ a' N _ n ). 

We can express H in terms of uncoupled fermion modes. First, given the translational 
invariance of the Hamiltonian, we introduce the Fourier modes 

1 ^ 2irk 

b k = -7^Yl e ~ i6kna n, ^ = V’ keZ - 

' n= 1 


Notice that by construction b *. = The Hamiltonian in 


H 


N-l 


£+~Y J , Ab- t ) 


k =0 



Q can now be written as 

AA 

T-J’ 


with 


which is real, 


L 

F k = Ae i6k \ 

l=—L 

L 

G k = B ie idkl , 

l=—L 


which is an odd function of k. i.e. G_ k = 

£ = 


—G k , and 


N- 

Y. f * 


k =0 


The matrix 



is Hermitian and satisfies 


( 2 ) 

(3) 

(4) 



Therefore we can introduce the dispersion relation A kl with A k+ N = A k , such that the 
diagonalization of R k is given by 


U k R k U\ = 



0 

—A- k 
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with Uk a unitary matrix. The order ambiguity in the eigenvalues is fixed if we demand 
that Af = (Afe + A_fc)/2 > 0. Thus the dispersion relation is univocally determined. 

We can now write the Hamiltonian in diagonal form by performing a Bogoliubov trans¬ 
formation to the new modes 



in terms of which H reads 

N -i , , \ 

H = £ + '£A t U k d t --\. 

k =0 ' ' 

A?; can be derived from the relations 

A A — Afc — A—k _ Fk — F-k _ J^A 
k ~ 2 “ 2 “ k ’ 

and 

A? = \j (U s ) 2 + |G t | 2 > 0, 

with 

vs — Fk + F^ k 
F “ = 2 ■ 

Finally we get 

An = ^(F/) 2 + |G*| 2 + Ff. (5) 

An eigenstate of H is determined by a subset of occupied modes K C {0,..., iV — 1}. 
If we denote by |0) the vacuum of the Fock space for the new modes, i.e. c4 |0) = 0, the 
stationary state 

|K) = n 4 |0> (6) 

fce K 

has energy 

z fceK z kg k 

In particular, the ground state will be obtained by filling the modes with negative energy 
(Dirac sea) 

ik) = n 4 io>, 

Afc<0 

N—l 

Fk = £ ~ d 

z k=0 


with energy 
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Note that if \F k \ < y/(F k ) 2 + \Gk\ 2 the dispersion relation is positive and the ground 
state of 77 is 10). On the contrary, if there are some momenta with \F k \ > yj (F k ) 2 + \Gk\ 2 
a Dirac sea develops, hence the ground state has occupied modes and differs from the 
Fock space vacuum. In this case, the ground state is not invariant under reflection in 
momentum space k -H- N — k. As we shall see in the next section, this fact has important 
consequences for the evaluation of the entanglement entropy and its behaviour near the 
boundaries of the critical region. 


III. ENTANGLEMENT ENTROPY AND CORRELATION MATRIX 


The main goal of this paper is to study the Renyi entanglement entropy for the ground 
state of a Hamiltonian of the kind analysed in the previous section. For the moment, 
however, we keep the general discussion by considering a general stationary state |K) as 
defined in (j6j) . 

Given a subsystem X of contiguous sites of the fermionic chain, with complementary 
set Y we have the corresponding factorization of the Hilbert space FL = Fix <E> FLy- If the 
system is in the pure state |K), the reduced density matrix for X is obtained by taking 
the partial trace with respect to FLy , px = Try |K) (K|. The Renyi entanglement entropy 
of X is defined by 

S«(X) = D-logins). (7) 

1 — a 

In order to compute S a (X) we take advantage of the fact that the state |K) satisfies 
the Wick decomposition property and, as it is well known HDum, the reduced density 
matrix can be obtained from the two point correlation function EH- 

For any pair of sites n, m G A", we introduce the following correlation matrix, 


(Vx) 


nm 


2 



-5 


nm 


/ 


2 (ai„a. 


m^n 


} 2 (a n a m ) 

^ ( a n a m ) ^ (aj,a m ) — 5 n 


Notice that Vx is a 2|X| x 2|X| matrix. 

Following mm one can show that the Renyi entanglement entropy reads 


S a {X) 


1 

2(1 - a) 


Tr log 


( I + Vx 


a 

+ 




a-\ 




( 8 ) 


where I is the 2|X| x 2|A| identity matrix. 
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From the numerical side, this formula supposes a drastic improvement of computational 
capabilities. In fact, the dimension of px is so the cost of calculating S a (X) grows, in 
principle, exponentially with the length of the interval. However, Vx has only dimension 
2\X\, allowing one to reach larger sizes of X. 

We can derive a useful way of implementing ([8]) by applying Cauchy’s residue theorem, 
so that 


S a (X) 


lim -— 

£ -s>o+ Am 


fa( 1 + £, A) 


dlog D x (X) 
dX 


dA, 


(9) 


where D x ( A) = det(AI — Vx), 


f a (x,y ) 


1 — a 


log 


y-v-y+y-v 


OL~ 


J 


( 10 ) 


and C is the contour depicted in Fig. [l] which surrounds the eigenvalues of Vx, all of 
them lying in the real interval [—1,1]. 



FIG. 1: Contour of integration, cuts and poles for the computation of S a (X). The cuts for the function 
f a extend to ±oo. 


In order to compute Vx we first rewrite it in the Fourier basis, 

1 JV " 1 
k =0 

with 

Or in terms of the Bogoliubov modes 




(4,cU) )u k -i. 
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We now compute the expectation value in the stationary state |K) associated to the 
occupation K C {0,1,..., N — 1}, so that 


di 


d ] 


(did-,)) = 


1 ~xA k ) 


o 


—k 


0 x K ( N ~ k ) 


(ii) 


where X K { k ) is the characteristic function of K. i.e. it is 1 or 0 according to whether k 
belongs to K or not. Now, introducing 


M k = Ul 


1 0 
0 -1 


u k = 


V(F k s y + \G k \* \ G k 


F k s G k 


-F k S 


we finally arrive at 


Gk = < 


( 12 ) 


-M k , if k E K and JV-fceK, 

-I, if keK and N - k g K, 

M k , if kg K and N -kg K, 

/, if k g K and N - k e K. 

Notice that the symbol of V\ at momentum k depends not only on the occupation number 
at k but also at N — k. This is a general fact that, to our knowledge, has been overlooked 
in the literature. We believe that there are two reasons for that. First, if F A (k) = 0 and 
we are interested in a ground state like that in ®E|, then we have K = 0. Hence any 
consideration on the occupation numbers is superfluous. The symbol is M k and that is 
all. The second situation corresponds to general stationary states, like in [4, [12j where 
the whole casuistics in the computation of the symbol would play a role. However in 
these papers, like in jT3] one has G k = 0, therefore M k is diagonal, the symbol can be 
considered scalar instead of a 2 by 2 matrix and the discussion above is not necessary. 

In this paper we are interested in computing S a (X) for the ground state. As it is well 
known, it contains information on the critical properties of the system. Following the 
previous discussion, it is determined by the occupation set 

K = {fc |A fc < 0}. 

Since A fc = Af + F' A , with Af = A^_ fc > 0 and Fjg = —Ffi_ k , applying the general 


expression in (12), we obtain 


Gk = < 


-A 

M k , 

I, 


if 

if 

if 


- ai > jy, 

- Af < F t A < Af, 


(13) 


Fl > Af. 
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Note that for the vacuum state, and due to the fact that A? > 0, one can not have that 


Afc and A ^-k are both negative implying that the first case in (12) never happens. 


In the so-called thermodynamic limit, N — y oo, l —> oo with Nji fixed, the previous 
IV-tuples, like A*,, Fk, Gk and the others, become 27r-periodic functions determined by the 
relation A (6k) = A Ok = 2nk/N and analogously for F(0), G(0 ), .... 

IV. BLOCK TOEPLITZ DETERMINANT AND ENTANGLEMENT 

ENTROPY 


In this section we compute an expression for D y(A) = det(AI — Vx) for large |A|. This 
enters in the formula ([9]) so that we can obtain the asymptotic behaviour of S a (X) for 
large \X\. 

There exists a generalization for block Toeplitz matrices of Szego theorem [TT] that 
allows to compute the linear dominant term of Dx( A). According to it, 


logD x (A)= A1 


logdet[A/ — G{0)\&9 + • • • , (14) 

where the dots denote sub-linear contributions. If we apply this general result to our case, 


the integrand with G(0) read from (13) becomes the simple result: 

c 



(A + l) 2 , 

if 

- A? > F k A , 


det[A I-Q{6)\ = < 

A 2 — 1, 

if 

~ Af < if < Af. 

(15) 


(A-l) 2 , 

if 

if > Af. 



After inserting this result into (14) and then (|9]) we immediately see that the linear, domi¬ 
nant contribution to the entropy vanishes. Therefore, in order to determine its asymptotic 
behaviour one is forced to compute the sub-dominant terms for the determinant that are 


hidden in the dots of (14). 


Case 1: F A {9 ) 2 < F s (9) 2 + \G(9 )\ 2 

We first discuss the case in which F A (0 ) 2 < F s (6) 2 + \G(0)\ 2 for all 0 6 (—7r, 7r]. The 
symbol G(6) is therefore continuous. 

This situation has been analysed before in the notable papers (5) [7] when F{6) = F(—9) 
and G(0) is imaginary. Their results can be straightforwardly generalised to our case 
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provided 

G(0) = -e^G(d). (16) 

The generalisation goes as follows. We first define the meromorphic complex functions 

L L 

= J2 A ‘ z ‘- s w = E Biz ‘• 

l=—L 1——L 

that are related to F and G, from Q and Q, by F{9) = d^e* 0 ) and G{6 ) = S(e* 0 ). We 
also introduce the symmetric and antisymmetric part of $, = ^($( 2 ) + $(^ -1 )), 

$ A (,z) = — <f>(z -1 )). Since A t = A_i and = —F>i, these functions satisfy 

$ s (^- 1 ) = W{z) = $ s {z), E (z- 1 ) = -E(z), 


and this, together with (16), implies 


3 (*) 



We now define the polynomial 

4 L 

P{z) = z 2L ($ 5 (E) + e^' /2 Z(z)) ($ s (z) - e^ 2 E(z)) = K Y[{z - Zj ), (17) 

3 =1 


with K = \(Al + Al) 2 — \Bl\ 2 . From the symmetries of <f> and 5 one can show that 
$' s (^) ±e^ ,//2 S(z) have real coefficients. Moreover, the roots of P(z) come in pairs Zj and 

A- 

We include the possibility of a global phase for completeness, but actually it is easy 
to get rid of it by performing a redefinition of the creation and annihilation operators. I 11 
fact, by introducing a! n = e"^ 4 an the new meromorphic functions are <&' = <3>, S' = e”^ 2 S, 
therefore 


(*), 


and the global phase is absent. 

We assume that P(z) is generic in the sense that it has 4L different single roots. Then 
the complex curve 

w 2 = P(z ) (18) 

defines a Riemann surface of genus g — 2L — 1. 
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As it is shown in [7], in the limit of large |X|, Dx(X) can be expressed in terms of 
analytic invariants of the Riemann surface C 9 . Namely, we make a choice of branch 
cuts in the z plane that do not cross the unit circle and a basis of fundamental cycles 
such that every a-cycle surrounds anticlockwise one of the cuts with the last L of them 
corresponding to the cuts outside the unit circle. Associated to this choice we have a 
canonically defined theta function, d : C 9 —> C, given by 


d(s) = d(s|n) = gfcrfl-nil+Wff n e Z g^ 

n 


where II is the g x g period matrix. Now, one has the following expression for the 
asymptotic value of the determinant (Proposition 1 of [7]): 


log D; y(A) 


\X\ log (A 2 - 1) + log 


&{P{\)e + |)$(/3(A)e 

0(lj* 


+ ... 


(19) 


where the dots represent terms that vanish in the large X limit. The argument of the 
theta function contains 

« A > = 2S l0 ^’ < 2 °) 

e G Z 9 which is a g dimensional vector whose first L — 1 entries are 0 and the last L are 1 
and the constant g-tuple r is determined once we fix which of the roots of P(z) are also 
roots of its first factor z L (<h s, (;j) + (see Ref. |7J for details). 


Case 2: F A {9) 2 > F s {6) 2 + \G{6)\ 2 

We move on now to discuss the other case, in which F A (9) 2 > F s {9) 2 + |G(d)| 2 for 
some open intervals in (—7r,7r]. The symbol now has discontinuities at the boundaries of 
the interval where F A {6) 2 = F s {6) 2 + \G(6)\ 2 . 

Up to now, in the presence of such discontinuities none of the known results on block 
Toeplitz matrices can be applied to obtain the asymptotic behaviour for the determinant. 
We here propose a procedure that covers this gap in the literature, so that one can compute 
the asymptotics of the determinant in this case. 

Before presenting the new procedure, let us recap some facts about Toeplitz matrices 
with scalar symbol with discontinuities and the corresponding Fisher-Hartwig theorem 
0 . 
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We therefore consider a genuine Toeplitz matrix, i.e. one with a scalar symbol g(0). 
Furthermore, we assume g{9) is a piecewise smooth symbol with jump discontinuities at 
9j, j — 1,..., J and lateral limits tj, fk at the discontinuities. Then from the discontinu¬ 
ities we get a logarithmic contribution to ZAv(A), whose asymptotic expansion reads 

logL> x (A) = ^y log[A — g(9)]d9 + log \X\ ^ ^log ( 21 ) 

where now the dots represent finite contributions in the large X limit, that can be 
computed explicitly but are not going to be necessary for us. The crucial fact is that 
contrary to the linear term (and also the constant one) the logarithmic contribution only 
depends on the behaviour of the symbol near its discontinuities, more concretely on its 
lateral limits. 

The previous observation can be used to compute the logarithmic term for a block 
Toeplitz matrix in some particular cases. Assume that a two-dimensional symbol Q{9) 
presents a discontinuity at say 9. We shall assume that the lateral limits, Q~ = 
lim e _ > .g-G(9) and Q + = lim . e _^§+G{9), commute. If the two matrices commute there is 
a basis where both are diagonal. Call and the corresponding eigenvalues at each 
side of the discontinuity. Inspired by the Fisher-Hartwig theorem for the scalar case, we 
argue that the contribution of this discontinuity to the logarithmic term of the determi¬ 
nant only depends on the value of Q{9) at each side of the jump or more concretely on 
its eigenvalues. Explicitly the contribution of the discontinuity to the coefficient of the 
logarithmic term of log Thy (A) is 


b = 


log 


A 


Ti 


— \ 2 


+ ^g 


A 


To 


— \ 2 


A — To 


( 22 ) 


47t 2 IV °A-Tj + 

Now if we have several discontinuities with commuting lateral limits there is a logarithmic 
contribution bj, like the one above, for each of them. Then, the asymptotic form of the 
block Toeplitz determinant up to finite contributions should be 


log D x ( A) 


m 

2vr 



j 

logdet[AJ — Q(9)]d9 + log |X| ^ bj H-. 

3 =1 


(23) 


This result, which is crucial for the computation of the asymptotic behaviour of the 
entanglement entropy for the general free fermionic model with finite range coupling, was 
not previously known in the literature. It is one of the main new contributions of this 
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work. It would be very interesting to extend this to the case when the two lateral limits 
do not commute. So far we have been unable to accomplish this goal. Fortunately, in 
order to study the fermionic chain that we consider in this paper, the stated result is 
enough. 

After this general discussion we proceed to study the asymptotic behaviour for the 
entanglement entropy in the ground state of the Hamiltonian ([Tj) when the symbol in the 
thermodynamic limit is not continuous. The latter can be written 


G(e) 


-/, if -A s (9) > F A (0), 

Af(0), if -A s (6) < F A (9) < A s (9), 

/, if F A (6) > A s (0). 


(24) 




The first thing one should notice is that if Q has a discontinuity at 6 it has also another 
one at —9. Therefore, discontinuities come in pairs except those at 6 = 0 or n. 

We will now analyse the type of discontinuities and their contribution to the logarithmic 
term. 

The first situation is when the two lateral limits are M{6) and ±/. Obviously, the two 
limits commute and we can therefore apply (22). The eigenvalues of M{6) are ±1 (indeed, 
det M{6) = —1 and Tr M{6) = 0), which implies that its contribution to the log |A"| term 


of (23) is 


h ^i = 4^2 ( lo S 


A + 1 
A — 1 


The second type of discontinuity is when the lateral limits are +/ and —in this case 
the two eigenvalues ±1 are different in both sides and we get a contribution 


bn = 2 b M i- 


Finally it is also possible that the matrix M{6) itself is discontinuous. This may 
happen when F s and G vanish for some values of 6 and at least one of them goes linearly 
to zero. In this case the two lateral limits have opposite sign and the contribution to the 
logarithmic term is 

bMM = 2 6 m /- 

From the above ingredients we can compute the asymptotic behaviour of the determi¬ 
nant for the four archetypical situations sketched in Fig. [2] 
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FIG. 2: Four archetypical discontinuities for the symbol G(9) of the correlation matrix. In a and b the 
dispersion relation A(8) is represented by the solid curve while the dashed curve depicts A(—6). In the 
plots c and d the solid curve stands for A s (8). The lines with arrows, right below the plots, mark the 
angle 9 where the discontinuities happen. 


In plot a we represent a double change of sign for the dispersion relation for positive 
values of 6 . We have four discontinuities of the kind MI , which one can deduce from 


(24). This is depicted with the lines with arrows below the plot. The contribution to the 
logarithmic term of log-Djy(A) is b a = AbMi- 

In plot b of Fig. [2] we consider the case in which A changes sign at 6 = 0 (the case 
6 = 7T is analogous). We now have two discontinuities of the type MI and one II with a 
contribution b b = 2 b mi + b n = 4 

In the plot c we represent the case in which A 5 vanishes at two symmetric values 
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of the angle, it produces two discontinuities of the type AIM and hence a contribution 
b c = 2b mm = 46 mi■ 

Finally it may happen that A s vanishes at 6 = 0 or 7 r in which case we have only a 
discontinuity of the type MM and a contribution 6 ^ = 26m/- Altogether we have 

6t = 2(2 n a + 27ib + 2 n c + rid)6 m/ = 6Vt bui- 


This leads to 


log Dx(X) = \X\ log(A 2 — 1) + log \X\br + • • • - 
Therefore, the asymptotic behaviour of the Renyi entanglement entropy is 

S a (X) — log |V| + • • ■ , 

where we have used the identity 

1 f , ,, , ,,db M i a + 1 
e->o+ 4m Jc x ’ ; dA 24a 


(25) 


(26) 


V. EXAMPLE I: DZYALOSHINSKI-MORIYA COUPLING 


As an application of the previous results we will study the XY spin chain with trans¬ 
verse magnetic held and Dzyaloshinski-Moriya (DM) coupling [9], [T7] . The Hamiltonian 
reads 

1 N 

Hbm = 2 + a)«+i + (t- A)Vn°n+i + s (<^+1 - <+!<) + ha z n ] . 

n= 1 


The coupling constants h, t, s and 7 are assumed to be real. 

As it is well know, a Jordan-Wigner transform allows to write this Hamiltonian in 
terms of fermionic operators. Namely 


N 


n= 1 


-66dm ^ ^ (6 + 7s ')+ (t 7s)(ll l Cl n —i T 7(^n^n+1 ^n^n+l) T 6 


Nh 


, (27) 


which can be described as a Kitaev chain with complex hopping couplings that break the 
reflection symmetry. It is also well known that the entanglement entropy of connected spin 
subchains coincide with its corresponding fermionic one and, therefore, for our purposes 
the two models are completely equivalent. 
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The Hamiltonian of (27) is a particular case of (|TJ) with L — 1 and then we can apply 
the last section’s results. The meromophic functions are in this case 


$ b (z) = tz + h+ tz~ 1 , $ A (z) = -is(z - z ), E(z) = ~/(z - z ). 

Therefore, the dispersion relation is 

A (0) = A s (6) + 2s sin 9, with A s (6) = yj(h + 21 cos 9) 2 + 4y 2 sin 2 9. (28) 

In the following we will fix the invariance under rescaling of the coupling constants by 
taking t — 1 . 


With a simple inspection of (28) we deduce that the system has no mass gap when 


A = s 2 — y 2 > 0, and (h/ 2 ) 2 — A < 1; (Region A), 

or when 

A < 0, and h = 2; (Region B). 

In Fig. [3] we depict regions A and B (actually a line) in the ( 7 , h) plane for a fixed s. 
In region A the dipersion relation becomes negative in some interval (Dirac sea) and 
therefore the energy is minimum when all these modes are occupied. It implies that the 
ground state corresponds to K 7 ^ 0 and we have discontinuities in the symbol Q(0). If 
h 7 ^ 2, it has four discontinuities of the type MI, which corresponds to the situation 
represented in Fig. [2] A. For h — 2 one of the zeros of the dispersion relation is at 9 = n 
and it corresponds to plot B in Fig. [2} Hence, we have two discontinuities of the type MI 


and one of the type II. In both cases, applying (26) we obtain an asymptotic behaviour 
for the entropy given by 

S a (X) = — 7 —— log |X| + ... 
ba 

In the points of region B A (9) is positive except at 6 = n where it vanishes. This 
implies that we have a single discontinuity of the type MM as we discussed before and it 
is represented in Fig. [2]D. The entropy in this region is 

s« ( v ) = °j^h° g |v| + ... 


In Fig. [f] we check numerically these scalings for the von Neumann entropy (a = 1). 
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h/2 



FIG. 3: Phase diagram for I?dm in the ( 7 , h ) plane for t = 1 and fixed s. The shaded Region A is 
gapless with central charge c = 1 (XX universality class) while the dashed line Region B has central 
charge c = 1/2 (Ising universality class). In the unshaded area, the Hamiltonian has a gap and there are 
conical curves with the same entropy. The dashed ellipses and hyperbolas depict some of them. 

We should discuss now the asymptotic behaviour of the entropy for the points in the 
( 7 , h) plane outside the critical regions A and B. I 11 this region the dispersion relation is 
always positive and 

G(0) = M(0 ), 

which besides is continuous, because A s (9) > 0. 

Notice that in our case 7 is real and therefore S(z) = S (z), hence we can follow the 
general procedure in terms of the Riemann surface and theta functions that we mentioned 
in the previous section. This has been actually carried out in [5j for the von Neumann 
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m 

FIG. 4: Numerical logarithmic contribution of |X| to the von Neumann entropy (a = 1) for two sets of 
couplings of Hdm- We have subtracted the constant term C\ = Si (100) — c/31og(100) with c = 1 or 
1/2. The □ correspond to s = 0.75, 7 = 0.5 and h = 0.5, so it belongs to Region A. The A represent 
s = 0.75, 7 = 1.5 and h = 2, so we are in Region B. The solid lines depict our prediction in each region. 


entropy and in [6] for the Renyi entropy. The final results can be stated in a very simple 
way. 

Let us introduce 

,, = 1 ~ (fa/ 2) 3 

7 2 

Now we should distinguish three regions outside the critical ones la, lb and 2. They are 
represented in Fig. [3j In every region the von Neumann entropy is constant with |X| and 
has a different expression. 


Region la: 0 < x < 1, y 2 > s 2 

1 — x 




log 


16 ^ 

Region lb: x > 1, y 2 > s 2 . 


+ 2(1+X) /(a/T^jT )I{^/~x 


7T 


+ log 2. 




log 


X 


+ 


2(1 +a; -1 ) 


7r 




x 


+ log 2. 
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• Region 2: x < 0, y 2 (l — x) > s 2 . 


S i = 


1 

12 



— x — x 



4(x — x ) i 

7r(2 — a; — 


VT 


x 



Here J(x) is the complete elliptic integral of the first kind 

_dy_ 

^(1 - y' 2 )(l - z 2 y 2 ) 

The hrst observation is that the entropy only depends on the constant x = (1 — 
(h/ 2 ) 2 )/ 7 2 . This fact was hrst noticed in [IS] where the case of s = 0 was discussed. 
The curves of constant x are ellipses in zones la and lb and hyperbolas in zone 2, they 
are ploted in Fig. [3] ft is clear that (if s = 0 and region A shrinks to a line) all curves 
intersect at h — 2, 7 = 0 which implies that the entropy does not have a well defined limit 
at that point. It has been called the essential critical point in [15], From the expression 
for the von Neumann entropy, one also observe a duality between the zones la and lb by 
transforming x into x _1 . In this way a curve in the region la is transformed into another 
one in lb while maintaining the same value for the entropy. One can show that the curves 
of constant entropy as well as the duality also hold for Renyi entropy. 

When we turn on the DM coupling (s 7 ^ 0) the critical line at 7 = 0 in Fig. [3] blows 
up and transforms into region A. The curves of constant entropy in the plane ( 7 , h) still 
persist, but only the portion of them outside region A. On the other hand the entropy at 
the essential critical point is well defined. As it is inside the critical region, it grows like 
the logarithm of the size of the interval and in the limit of large |X| it becomes 

q r l f 4s 2 A 

s “ {x) = ~&r log (r+A' Y| ) +It " 



where X a is a constant independent of s whose expression is 


T =_ 

-L'OL 


m 


d/q(l, A) 
dA 


log 


r (l/2-ff(A)) 

L r(l /2 + /3(A) 


dA, 


where F represents the Euler gamma function, f a is defined in (10) and /3(A) in (20). 


Another special point is 7 = 1, h — 2 where, as was noted in [9] by Kadar and 
Zimboras, //dm is Kramers-Wannier self-dual. In that case, the correlation matrix reduces 
to a Toeplitz matrix with scalar symbol and employing the Fisher-Hartwig theorem they 
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can obtain not only the logarithmic contribution but also the constant term. For the von 
Neumann entropy they have 

o J r | 1 °g( 2 l A 1) + ^ lo g( 1 - s_2 )+ :z: i; S>1 

|pog(4|X|) + f; »<1. 

Note that the logarithmic term agrees with our general result. 

Another interesting property of this theory is that when we approach the critical zone B 
from la or 2, the entropy saturates at a larger and larger value that grows logarithmically 
with the correlation length (or the inverse of the mass gap). This is the expected behaviour 
that has been predicted from the properties of conformal field theories [TB] and has been 
verified numerical and analytically in many models. 

As for the boundary of region A, it can be reached from inside the region or from 
outside. These two limits are completely different. Of course, inside the critical region 
the entropy diverges logarithmically with the length of the interval and, what matters 
for us, the coefficient of this term is constant throughout all the region. Contrary to this 
behaviour, the constant coefficient in the asymptotic expansion of the entropy does change 
inside region A and it indeed diverges with negative values when we tend to any boundary. 
If, we now study the case in which we approach zone A from any of the non critical ones 
(la, lb or 2) we find that the entanglement entropy of the interval saturates in the limit 
at a finite value, which is different at any point of the boundary but independent of the 
path (always inside the non critical region) that we follow to reach the point. Note that 
this behaviour is anomalous in the sense that, as we mentioned before, one expects an 
exponential divergence of the entropy with the inverse of the mass gap. Recall that in 
the non critical regions la, lb and 2 the ground state is the Fock space vacuum for the 
Bogoliubov modes; the same occurs in the critical region B. On the contrary, in region 
A the ground state has all Bogoliubov modes with negative energy occupied (the Dirac 
sea). The above behaviour of the entropy near the transition could be a sign of this 
discontinuity in the ground state. In fact, we find a similar issue in the XX spin chain 
(i.e. when s,y = 0). In that case, the entropy always vanishes when h > 2, it does 
not diverge with the inverse of the mass gap when we approach the critical point h — 2, 
and the ground state is the Fock space vacuum. On the contrary, when h < 2, the von 
Neumann entropy scales like 1/3 log |X| while the constant correction diverges to negative 
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values when we approach the essential critical point. In this region a Dirac sea develops 
and the ground state has occupied all the modes with negative energy. The situation is, 
therefore, very similar to what we obtained for the DM model. 

From a global point of view, this behaviour in the border of region A can be interpreted 
as a blow up of the essential critical point of Ref. SSI- The different possible values for 
the limit of the entropy at that point (for s = 0 ) are obtained at different points of the 
boundary of region A for s 7 ^ 0 and the singularity at h = 2 , 7 = 0 disappears in this 
case. 

VI. EXAMPLE II: KITAEV CHAIN WITH LONG RANGE PAIRING 

As a further application of the previous results we shall discuss the case in which the 
couplings extend throughout the whole chain. The example is adapted from [TSJ and it 
serves to illustrate how to use the tools of the preceding sections to determine the scaling 
behaviour of the entropy in a theory with long range couplings. 

The Hamiltonian represents a Kitaev chain with power like decaying pairing, i. e. 

N N N/2 

Hk ^ ^ T & n +\Cl n T H CL^CLf^j T ^ ^ ^ l\l\ ^ (®n®n-|-£ ^n^n+l) ^ t (29) 

n =1 n =1 l=—N/2 

where the exponent ( > 0 characterises the dumping of the coupling with distance. Its 
value will happen to be critical to determine the scaling behaviour of the entropy. 

Following the discussion in section 4 and taking the thermodynamic limit, N oo, 
we define 


$ 5 (^) = $(z) = z + h + z 1 (30) 

OO 

3c(*) = - z ~ l ) l ~ ( = Li t0) - Li c0 _1 ); ( 31 ) 

i=i 

where Li^ stands for the polylogarithm of order (. This is a multivalued function, analytic 
outside the real interval [ 1 , oo) and, what is going to be most important for us, has a finite 
limit at z = 1 for ( > 1 while it diverges at that point for £ < 1 . 

If we introduce now G^(6) = S^e* 0 ) which is a purely imaginary function, the disper¬ 


sion relation reads 


A(0) = y(h + 2cos0) 2 + |G C (0)| 2 , 
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and vanishes at 9 = n for h = 2 and at 9 = 0 for h = —2 and ( > 1. These are the two 
instances in which the mass gap is zero. An interesting feature of this model is that, as 
we shall see, it may have logarithmic scaling of the entropy even when it has non zero 
mass gap. 

Notice that G(9) = —G(8) and therefore the formalism of section 4 can be applied, we 
also have A A {6) = 0 which corresponds to cases c or d in section 4. Hence, to uncover the 
possible logarithmic scaling of the entropy we must simply look for discontinuities of the 
matrix 

1 ( h + 2cos9 Gr(9 ) 

M (9) = -tTq\ ( 

H d ) ^ -G c (9) —h — 2 cos6* 

One source of discontinuities comes from the zeros of the dispersion relation, this 
corresponds to the cases c or d discussed in section 4. In our model the only zeros appear 
at 9 = 0 for h = —2 or 9 — 7T for h — 2 and the two lateral limits of M at the discontinuity 
are ±cr y . Therefore they contribute to the effective central charge for the scaling of the 
entanglement entropy with c = 1/2. Note that this contribution has its origin in the 
absence of mass gap, i. e. it is connected to the universality class of a conformal held 
theory. 

The other possible discontinuities are related to the divergences of G{9) at 9 = 0; 
these happen for ( < 1 and any value of h. The two lateral limits are also ±a y and 
the discontinuity contributes with c = 1/2 to the scaling of the entropy. If h — 2 this 
contribution must be added to the one coming from the zero of the dispersion relation at 
9 = 7r, therefore the total central charge is c = 1. When h = —2 such an addition does 
not happen as the two discontinuitites are actually the same and the total central charge 
is c = 1/2. 

All these results are summarized in Fig. [5] where the different regions, according to 
the value of c, are shown. The white region of central charge c = 0 corresponds to ( > 1 
and h ^ ±2; in black we represent the region with c = 1 i. e. h = 2 and £ < 1; the rest 
in gray corresponds to c = 1/2. 

Our analytic results are very much compatible with the numeric hirelings of [18J, where 
the regions we have determined appear smeared somehow. We interpret this fact as the 
consequence of finite size corrections to the thermodynamic limit that we are presently 
considering. Note that for comparison with Ref. pT8j one should bear in mind that their 
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o 

FIG. 5: Plot of the different regions in the (h, C) plane according the the coefficient of the logarithmic 
scaling for the entanglement entropy. White area stands for c = 0, grey region for c = 1/2 and the black 
one for c = 1. 


-4 -2 0 2 4 
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parameters a and fi correspond to ( and h/2 in our paper. 

Another interesting remark is that the results above can be applied to any family of 
long range couplings with the same asymptotic behaviour for large l and constant sign. 
In fact, our discussion was based on the limiting behaviour of G at 6 = 0 and the latter 
is governed solely by the asymptotic behaviour of the power series. 


VII. CONCLUSIONS 

In this paper, we extend further the understanding of entanglement entropy in spin¬ 
less fermionic chains considering translational invariant Hamiltonians with finite range 
couplings breaking the fermionic number, the reflection and the charge conjugation sym¬ 
metries. The first violation implies to deal with block Toeplitz determinants while the 
second one can create Bogoliubov modes with negative energy. The latter breaking leads 
to discontinuities in the symbol of our block Toeplitz matrix. To our knowledge, this pos¬ 
sibility has not been considered in the literature. Here, we have performed a systematic 
analysis of their physical origin, nature and contribution to the entanglement entropy of 
a set of contiguous sites in the thermodynamic limit. For the last point, we have carried 
out an heuristic study based on the fact that the contribution of each discontinuity only 
depends on the value of the symbol at each side of the jump. Thereby, we conclude that 
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each of them adds a logarithmic term with the size of the interval to the entropy. This 
method also allows to determine its coefficient which informs about the universality class 
of the critical theory. From a mathematical point of view, our result requires that the 
two lateral limits of the symbol at the discontinuity commute. Although this is enough 
for our physical problem, it will be a very interesting mathematical question to study the 
case in which the two lateral limits do not commute. 

When the reflection symmetry breaking does not produce Bogoliubov modes with 
negative energy, the expressions found by Its, Mezzadri and Mo can still be applied. Here 
we have also shown that it can be extended to certain cases where the charge conjugation 
symmetry is broken. 

We have applied this general analysis to a XY spin chain with transverse magnetic 
field and Dzyaloshinski-Moriya coupling. This system has two critical regions as was 
previously pointed out in [9] and m Here we determine that one region belongs to the 
Ising universality class while another one region corresponds to the XX universality class. 
The latter can be seen as a blowing up, produced by the DM coupling, of the XX critical 
line of the XY spin chain with only magnetic coupling. 

We have also tested our formula in a theory with longe range couplings, a Kitaev 
fermionic chain with power like decaying pairing. In this case, the logarithmic scaling of 
the entropy is affected by the dumping of the coupling with the distance. For sufficiently 
small values, it produces a discontinuity in the symbol of the correlation matrix. As 
a result, for zero mass gap and positive chemical potential, the univesality class of the 
critical theory changes from the XX class (for small values of the dumping) to the Ising 
class (when the dumping is sufficiently large). In addition, this discontinuity implies that 
the entropy may have a logarithmic scaling when the mass gap is not null. Our analytical 
study is in a great deal agree with the numeric results obtained in [TH]. The difference is 
that the regions we have precisely delimitated here, appear blurred in [IB]. We interpret 
this fact as a consequence of finite size corrections to the thermodynamic limit we consider 
here. 
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